#ifndef __ROUTINE_LINEAR_RESPONSE__
#define __ROUTINE_LINEAR_RESPONSE__

#include "stdio.h"
#include "stdlib.h"
#include "complex.h"
#include "math.h"

#include "string.h"

#define c_vacc (137.036)       // in a.u. units
#define PI (3.1415926535897932)

/**
 * 1. load the ground state of molecule from data file
 * 
 */

void load_ground_state(char* file_name, int cutstate, int st_num,     // infos
		       double *e_values, double *occ_num,  // orbital informations
		       double *rho_0, double *wfs,         // density and orbitals
		       double *v_0_ext, double *v_0_H, double *v_0_xc,
		       int GsNx, int GsNy, int GsNz,
		       int Nx, int Ny, int Nz);


/**
 * 2. evaluation of the value of rho and j with new P
 * 
 */

void eval_rho_j(double complex *del_rho, 
		double complex *del_j_1, 
		double complex *del_j_2, 
		double complex *del_j_3, // return values
		double complex *P, double complex theta, 
		double *wfs, double *gradient_wfs, 
		double *e_values, double *occ_num, 
		int N_sample, int occ, int unocc, 
		int Nx, int Ny, int Nz);

/**
 * update the right hand side of Casida's equation
 * 
 */

void update_rhs(double complex *F,             // output
		double complex *AEM_1, 
		double complex *AEM_2, 
		double complex *AEM_3, 
		double complex omega,
		double *wfs, double *gradient_wfs, 
		int N_sample, int occ, int unocc, 
		int Nx, int Ny, int Nz,
		double hx, double hy, double hz,
		double *Intl_coef, double *QD_emass);

/**
 * update the matrix of Casida's equation 
 * where S*p := S1*p + S0.
 * It is much cheaper than sparse matrix multiplication!
 */

void update_matrix(double complex *S0, double complex *S1,
		   double complex *AEM_1,  double complex *AEM_2, 
		   double complex *AEM_3,           // EM field
		   double complex *del_rho, double complex *del_j_1, 
		   double complex *del_j_2, double complex *del_j_3,   
		   double complex omega,  double *rho_0,
		   double *drho_0_x,double *drho_0_y,double *drho_0_z,
		   double *wfs, double *gradient_wfs,
		   double *e_values, double *occ_num,  
		   double complex *fxc_kernel, double complex *fxcL0, 
		   double complex *fxcL, double complex *fxcT,
		   double complex *eta_xc, double complex *xi_xc,
		   double complex *gradient_eta_xc, 
		   double complex *hessian_eta_xc,     // help functionals 
		   int N_sample, int occ, int unocc,  //2.Quantum config 
		   int Nx, int Ny, int Nz,         // 3.computational mesh
		   double hx, double hy, double hz,
		   double *Intl_coef,
		   double *QD_emass, double *QD_epsilon);


int BiCGStab_Casida(double complex *S0, double complex *S1,
		    double complex *x, double complex *b,
		    int Msize, int precond, int MAXIT, double TOL);

#endif
